* Program to construct bootstrap percentiles of noise-signal ratio from real-time gap
* data over the 1980-1992 sample period.
* Note:  blockwidth and nobs are set by hand, and this program requires the number
* of blocks to divide evenly into the sample for the circular MBB procedure.


* This program:  Uses a block width of 4 and the RMSE-based noise-signal ratio


calendar 1947 1 4
allocate 2013:04+5000

seed 4892

open data table1_boot_data.rat
data(format=rats) / gb_gap_rt gb_gap_gdp94q3

* GB_GAP_RT:       Real-time gap estimate
* GB_GAP_GDP94Q3:  Final gap for 1980-1992 period


* Put the relevant data series into a vector (note that
* the make command will actually yield an array...)

make rtgap_vec 1980:01 1992:04
# gb_gap_rt
disp %rows(rtgap_vec); * Check that this equals NOBS, below

make fingap_vec 1980:01 1992:04
# gb_gap_gdp94q3
disp %rows(fingap_vec); * Check that this equals NOBS, below

compute gaprev = fingap_vec-rtgap_vec
compute sqgaprev = gaprev.^2.0

compute blockwidth=4
compute nobs=52
compute numblks=13
compute numdraws=5000

declare vector xxx(numblks)
declare vector[integer] entries(nobs)
declare vector b_sqrevs(nobs)
declare vector b_fin(nobs)
declare vector n_s_ratios(numdraws)



do drawiter=1,numdraws


ewise xxx(i)=%raninteger(1,nobs)

* Use these two lines of code (with NUMDRAWS=1) to test the calculations (they just return the original series)
* compute qqq=%seqa(1,blockwidth,nobs)
* ewise xxx(i)=qqq(i)


do iter=1,numblks
  do subiter=1,blockwidth
    compute entries(((iter-1)*blockwidth)+subiter)=%clock(fix(xxx(iter))+(subiter-1),nobs)
  end do subiter
end do iter

do iter=1,nobs
  compute b_sqrevs(iter)=sqgaprev(entries(iter),1)
  compute b_fin(iter)=fingap_vec(entries(iter),1)
end do iter

compute rmse=sqrt(%avg(b_sqrevs))
compute sqdev=(b_fin-%avg(b_fin)).^2.0
compute nnn=%rows(sqdev)
compute nscale=float(nnn)/float(nnn-1)
compute stddev=sqrt(nscale*%avg(sqdev))

compute n_s_ratios(drawiter)=rmse/stddev

end do drawiter



set nsr_series 1 numdraws = n_s_ratios(t)
statistics(fract) nsr_series

